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ABSTRACT 

Spatially resolved studies of star forming regions show that the assumption of spherical 
geometry is not realistic in most cases, with a major complication posed by the gas being 
ionised by multiple non-centrally located stars or star clusters. Geometrical effects including 
the spatial configuration of ionising sources affect the temperature and ionisation structure 
of these regions. We try to isolate the effects of multiple non-centrally located stars, via the 
construction of 3D photoionisation models using the 3D Monte Carlo photoionisation code 
MOCASSIN with very simple gas density distributions, but various spatial configurations for 
the ionisation sources. Our first aim is to study the resulting temperature structure of the 
gas and investigate the behaviour of temperature fluctuations within the ionised region. We 
show that geometry affects the temperature structures in our models differently according to 
metallicity. For the geometries and stellar populations considered in our study, at intermedi- 
ate and high metallicities, models with ionising sources distributed in the full volume, whose 
Stromgren spheres rarely overlap, show smaller temperature fluctuation than their central ion- 
isation counterparts, with fully overlapping concentric Stromgren spheres. The reverse is true 
at low metallicities. Finally the true temperature fluctuations due to the stellar distribution (as 
opposed to the large-scale temperature gradients due to other gas properties) are small in all 
cases and not a significant cause of error in metallicity studies. 

Emission line spectra from H II regions are often used to study the metallicity of star- 
forming regions, as well as providing a constraint for temperatures and luminosities of the 
ionising sources. Empirical metallicity diagnostics must often be calibrated with the aid of 
photoionisation models. However, most studies so far have been carried out by assuming 
spherical or plane-parallel geometries, with major limitations on allowed gas and dust den- 
sity distributions and with the spatial distribution of multiple, non-centrally located ionising 
sources not being accounted for. We compare integrated emission line spectra from our models 
and quantify any systematic errors caused by the simplifying assumption of a single, central 
location for all ionising sources. We find that the dependence of the metallicity indicators on 
the ionisation parameter causes a clear bias, due to the fact that models with a fully distributed 
configuration of stars always display lower ionisation parameters than their fully concentrated 
counterparts. The errors found imply that the geometrical distribution of ionisation sources 
may partly account for the large scatter in metallicities derived using model-calibrated empir- 
ical methods. 
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1 INTRODUCTION 

The ability to measure accurate chemical abundances in H II re- 
gions in our own and other galaxies is vital for our understanding of 
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their chemical evolution. Emission lines emitted by the nebular gas 
photoionised by massive stars provide us with powerful metallicity 
indicators for near and intermediate redshift galaxies. The major 
complication is posed by the critical task of determining the physi- 
cal parameters of the nebula, in particular the electron temperature, 
T e . Failure to achieve realistic estimates may lead to gross errors in 
the final abundance determination. Temperature fluctuations within 
the nebula may be an important cause of error. Collisionally excited 
lines (CELs), which are routinely used in abundance studies, are 
naturally weighted towards hotter regions and may therefore lead 
to underestimating of the real abundances. Recombination lines 
(RLs) are less affected by errors in temperature determinations, and 
should in theory yield more accurate results. One of the outstanding 
problems in nebular astrophysics is the discrepancy between abun- 
dances and electron temperature estimates obtained from CELs or 
RLs in Planetary Nebulae (PNe) (e.g. Rubin et al. 2002; Wesson, 
Liu & Barlow 2005; Ercolano et al. 2005; Liu et al. 2006; Peim- 
bert & Peimbert 2006 and references therein), galactic and Magel- 
lanic H II regions (e.g. Peimbert, Peimbert & Ruiz 2000; Peimbert 
2003; Tsamis et al. 2003; Esteban et al. 2004; Garcfa-Rojas et al. 
2005, 2006) and extragalactic H II regions (e.g. Peimbert & Pe- 
imbert 2003; Peimbert, Peimbert & Ruiz 2005). The cause of this 
discrepancy is still uncertain; temperature fluctuations and chemi- 
cal inhomogeneities in the gas have both been advocated to explain 
the discrepancy. Photoionisation models including chemical inho- 
mogeneities have been successful in matching the observed CEL 
and RL spectra of some PNe (e.g. Ercolano et al. 2003b) and H II 
regions (Tsamis & Pequignot 2005). However the discrepancy be- 
tween the abundances derived from recombination lines (RL) and 
those derived from CELs has been also been attributed to temper- 
ature fluctuations not caused by abundance inhomogeneities (e.g. 
Garcfa-Rojas et al. 2004). The t 2 parameter (Peimbert, 1967; Peim- 
bert & Costero 1969; Peimbert 1971), introduced in order to quan- 
tify such fluctuations, can be derived empirically by the compari- 
son of temperatures obtained from the depth of the B aimer Jump 
to those obtained from CELs. Although the situation for H II re- 
gions is notoriously less worrying than for PNe, the problem still 
remains that empirically determined t 2 values are larger than those 
indicated by chemically homogeneous photoionisation models, in- 
dicating that the causes of the proposed temperature fluctuations 
are still not understood. 

Aside from temperature fluctuations the task of obtaining reli- 
able abundance estimates from CEL spectra of H II regions is fur- 
ther complicated by the fact that often the complete set of emission 
lines needed for a direct measurement of the electron temperatures 
is not available for all sources; the 4363 A line of [O III] and the 
5755 A line of [N II] are weak and therefore not detected if the spec- 
tra do not have good signal-to-noise or if the metallicities are high, 
implying low electron temperatures. A large effort has been made 
to provide metallicity and temperature indicators based on various 
combinations of strong lines (e.g. Pagel et al. 1979; Alloin et al. 
1979; Mc Gaugh 1991; Storchi-Bergmann et al. 1994; Vflchez & 
Esteban 1996; Van Zee et al. 1998; Pilyugin 2001a; Pettini & Pagel 
2004; Perez-Montero & Diaz 2005). Hiigele et al. (2006) have re- 
cently proposed a new methodology for the production and calibra- 
tion of empirical relations between the different line temperatures 
based on observational data alone. However the data sets available 
are still too limited to provide usable indicators. Currently, calibra- 
tions of metallicity indicators and ionic temperature relations are 
generally not based on observational data alone, but make use of 
grids of photoionisation models, run for a range of metallicities and 
ionisation parameters (e.g. McGaugh 1991, Chariot & Longhetti 



Table 1. Characterising model parameters. Lbol is the total bolometric lu- 
minosity of all stars included. R s is the Strbmgren radius in the case when 
all stars are located at the centre of the sphere/shell. The total number of ion- 
ising photons and the hydrogen number density are constant for all models 
and are, respectively, Q H o = 3.80-10 50 sec -1 and Nn = 100 cm' 3 . All 
shell models have an inner radius of 2.8- 10 19 cm. 
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Figure 1. The left panel shows a 3D representation of the Strbmgren sphere 
distribution for case F, plotted as the iso-surfaces where the ionisation frac- 
tion of hydrogen is 0.95. The adjacent right panel shows an average projec- 
tion map of the ionic abundance of H+ . 



2001, Kewley & Dopita 2002). Given the vast parameter space gen- 
erally under investigation, most studies so far have been carried out 
with spherically symmetric or plane parallel geometries, with major 
limitations on the allowed density distribution and with the spatial 
distribution of the ionising sources not being investigated. Nearby 
H II regions, however, show complex structures in the distribution 
of gas and stars, which are often intermixed. 

In this paper we carry out a theoretical investigation of the 
importance of the effects due to the spatial distribution of the ionis- 
ing sources via the construction of a number of 3D photoionisation 
models, using the MOCASSIN code (Ercolano et al. 2003a, 2005) 
for simple gas density distributions and three spatial configurations 
for the ionisation sources. 

Our modelling strategy is described in Section 2. Temperature 
fluctuations are estimated for our models by computing theoreti- 
cal t 2 and To values (Peimbert, 1967) for each of our models and 
the results are given and discussed in Section 3. In this Section we 
also compare integrated emission line spectra from such configu- 
rations to search for systematic errors which may be caused by the 
simplifying assumption of a single, central location for all ionis- 
ing sources. We also test the robustness of a number of commonly 
used ionic temperature relations. A summary of our main results is 
presented in Section 4. 
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Table 2. Model grid legend. H number density, Nfj = 100 [cm 3 ] for all models. Figures in this paper rely on colour, please refer to on-line version. For 
black and white copies black is black, blue is dark-gray, red is medium gray and green is light gray. 
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2 MODELLING STRATEGY 

The aim of the current investigation is to uncover possible sys- 
tematic differences on the temperature structure and emission line 
spectra of nebulae ionised by a centrally concentrated set of stars 
or clusters versus those ionised by the same set of sources which 
are randomly distributed within the half or full volume, whose re- 
spective Stromgren spheres do or do not overlap. In particular we 
want to test which of the stellar configurations considered produces 
the largest temperature fluctuations and of what magnitude. At this 
stage we do not attempt to assess the absolute strengths and weak- 
nesses of one metallicity indicator or ionic temperature relation 
over others, a task that may only be carried out by a systematic 
investigation of the large parameter space. In an attempt to iso- 
late the effects of the stellar distribution from those due to the gas 
density distribution, we consider two extremely simple geometries 
- a homogeneous spherical volume and a homogeneous spherical 
shell, both of constant hydrogen number density, Nn = 100 cm -3 . 
The shell models have inner radii of 2.8- 10 19 cm and Stromgren 
radii (corresponding to the case when all stars are located at the 
centre) as listed in Table[T] The total number of ionising photons is 
constant for all models and is Q H a = 3.80-10 50 sec -1 . 

For each density distribution we are interested in a compari- 
son between the centrally concentrated and the half and fully dis- 
tributed source cases. In the remainder of this paper we will refer 
to models ionised by a central concentration of stars as models C, 
those ionised by the same set of stars distributed in the half and full 
volumes will be referred to as models H and F, respectively. Case C 
is equivalent to models ionised by a single central source with a to- 
tal bolometric luminosity equal to the sum of the bolometric lumi- 
nosities of the individual sources (as given in Table[TJ and a spectral 
shape given by the superposition of the spectra of all sources. These 
models could also be performed using a ID code, with a spherically 
symmetric gas density distribution. Similarly a ID code could also 
be employed for case F models where the Stromgren spheres do 
not generally overlap^; the left panel of Figure Q] shows a 3D rep- 

1 Due to the stochastic nature of the ionisation source distribution, it is pos- 



resentation of the Stromgren sphere distribution for case F, plotted 
as iso-surfaces where the ionisation fraction of hydrogen is 0.95. 
The adjacent panel shows a projected map of the ionic abundance 
of H + . In the intermediate case, case H, the Stromgren spheres of 
most H II regions partially overlap; the determination of the radia- 
tion field in the overlap region is not possible without the applica- 
tion of a 3D code. As we are interested in a comparison of the three 
cases (C, H and F), self-consistency is crucial and for this reason 
we have run all models with the same 3D code, MOCASSIN. 

Our models include five metallicities (Z/Zq = 0.05, 0.2, 0.4, 
1.0, 2.0). The "solar" abundance set uses the values from Grevesse 
& Sauval (1998) with the exception of C, N, O abundances which 
are taken from Allende Prieto et al 2002, Holweger 2001 and Al- 
lende Prieto et al 2001, respectively. These were scaled to lower 
and higher metallicities considering the empirical abundance trends 
observed in H II regions by Izotov et al. (2006). In order to max- 
imise the effects of spatially distributed sources whose Stromgren 
spheres may overlap (completely or partially), or be totally inde- 
pendent, we consider the limiting case of an ionising set composed 
of two stellar populations, a 37 Mq and a 56 M© population, with 
half of the total ionising photon output [s -1 ] being emitted by each 
population. These two stellar masses were chosen as they have very 
different Qhc°/Qh° ratios, where Q H o is the total output [s -1 ] of 
H-ionising photons (energy > 1 ryd) and Qh c ° is the total output 
[s ] of He-ionising photons (energy > 1.8 ryd), and are likely to 
produce the largest effects on the temperature structure and sharp- 
ness of the ionisation front. The ionising spectra for single-mass 
stars were computed with the STARBURST99 spectral synthesis 
code Leitherer et al. (1999) with the up-to-date non-LTE stellar at- 
mospheres implemented by Smith, Norris & Crowther (2002), us- 
ing single isochrones for the appropriate stellar masses. The models 
were calculated at metallicities consistent with the nebular gas and 
were obtained for an instantaneous burst, at an age of 1 Myr. At so- 
lar metallicities for the 37 Mq stars, the stellar atmosphere models 



sible that a small fraction of them may actually have partially overlapping 
Stromgren spheres. 
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emit 32.3% of luminosity in the H-ionising continuum and 7.7% 
in the He-ionising continuum, while for the 56 Mq stars these per- 
centages are 47.9% and 13.7%, respectively. The exact percentages 
vary with stellar metallicity, nevertheless the values above are given 
as a guide to appreciate the different spectral hardness of the two 
populations. 

Some other defining parameters of our models, together with 
their nomenclature and associated symbols are given in Table[2] 



2.1 The3Dphotoionisationcode:MOCASSiN 

The 3D photoionisation code MOCASSIN (Ercolano et al. 2003a, 
2005) uses a Monte Carlo approach to the radiation transport prob- 
lem and it is therefore completely independent of geometry and 
density distribution. Both the stellar and diffuse components of the 
radiation field are treated self-consistently, without the need of ap- 
proximations. Multiple ionisation sources can be located at arbi- 
trary positions in the simulation grid with the only limit being im- 
posed by computing resources. The atomic data used is frequently 
updated and include sets of energy levels, collision strengths and 
transition probabilities from Version 5 of the Chianti database 
(Landi et al., 2005) and the improved H I, He I and He II free-bound 
continuous emission data recently published by Ercolano & Storey 
(2006). A public version of the FORTRAN 90 code, which is fully 
parallelised using the Message Passing Interface (MPI) libraries, 
can be obtained from B. Ercolano. Version 2.02.38 was used for the 
models presented in this work. We simulate one quadrant of each 
model, using the technique described and tested by Ercolano et al. 
(2003b), whereby the positive x-y, y-z and x-z planes intersecting 
the z-, x- and y- axes respectively at zero, act as mirrors, reflect- 
ing the incoming photons back into the simulated cube. The full 
volume is finally described by 10 6 cells for the shell models and 
by 125000 cells for the spherical models ionised by 240 sources. 
The number of energy packets used varies during the course of 
our simulations, but typically 1-10 million packets are sufficient for 
our grids to reach convergence within 10-20 iterations. We experi- 
mented with higher resolution grids and a larger number of energy 
packets and found our results to be virtually unaffected. 



2.2 Validation of our models 

As stated above, the aim of this paper is not to provide new cali- 
brations to abundance diagnostic ratios nor to assess their absolute 
accuracy. It is therefore not in our intention to create models that 
fit any particular observations. While trying to maximise the effects 
of a complex ionising field distribution, it is however still necessary 
to ensure that the ionisation and temperature structures and, hence, 
emission line ratios we obtain from our models are in the range of 
those observed in nature. In Figure|2]we plot our results in a number 
of line ratio diagrams, including those following Veilleux & Oster- 
brock (1987) and Osterbrock, Tran & Veilleux (1992), that show 
the variation in [N II] and [S II] excitation parametrised in terms 
of [O m]/H/3 (top panels) and those following Mc Call, Rybski & 
Shields ( 35) showing the variation in terms of R230- Similarly, in 
Figure[3]we plot the H II region ionisation sequence parametrised in 
terms of oxygen and sulphur line ratios 0. The grey dots represent 
giant H II regions in spiral galaxies, taken from Garnett & Kenni- 
cutt (1994), Garnett et al. (1997), van Zee et al. (1998), Bresolin 



R23 = ([O III] 5007,4959 + [O II] 3726,3729)/H/3 

r)' = ([O n] 3726,3729/[0 III] 5007,4959)/([S II] 6716,673 1/[S III] 9069,9532) 



et al. (1999, 2004, 2005), and Kennicutt et al. (2003), spanning 
metallicities from 12 + log(0/H) ~ 8 to ~ 8.8. Figure |2] and the 
top right panel of figure [3] also include a few points from metal- 
poor emission line galaxies from Izotov et al. (2006), which extend 
to 12 + log(CVH) ~ 7.2; it was not possible to use this data for 
the other three panels of figure [3] due to the lack of simultaneous 
detection of all the oxygen and sulphur lines needed. In general, 
our models fall within or near the locus of observed H II regions, 
with the exception of our lowest metallicity (Z/Zq = 0.05) and 
highest metallicity (Z/Zq — 2) models, which is not surprising 
given that the observational sets available did not include many data 
points at such extreme values of Z; furthermore we did not attempt 
any systematic variation of the ionisation parameter with Z, while 
observations of giant H II regions argue that the ionisation param- 
eter decreases as Z increases (see e.g. Mc Gaugh 1994). However 
this does not affect the achievement of our aims. 



3 RESULTS 

The temperature structure and emission line spectra obtained by our 
models are analysed here in detail. Temperature fluctuations, which 
may introduce errors in the empirical calculations of abundances, 
are examined. The robustness of commonly used abundance diag- 
nostics and ionic temperature relations are also tested against pos- 
sible errors introduced by geometrical effects. 



3.1 Temperature structure 

The volume integrated electron temperatures weighted by a num- 
ber of commonly observed ionic species are plotted in Figure|4]for 
metallicities of Z/Zq = 0.05, 0.2, 0.4, 1 and 2. The values were 
obtained according to 



f N c N(X +[ )T c dV 
fN c N(X+ l )dV 



(1) 



The models are represented by circles and squares of various 
colours and sizes as described in Table[2] 

At a first glance some trends are apparent. At solar and high 
metallicities (i.e. Z/Zq > 1) the electron temperatures in 'high' 
ionic species zones (e.g. O III, S III, Ar III and N III) are higher 
for models in cases H and F. The opposite is true for electron tem- 
peratures in 'low' ionic species zones, such as O II and N II. This 
implies a shallower temperature gradient for case F and H when 
compared to case C, since at these metallicities the highest nebu- 
lar temperatures are reached in the low ionisation zones which are 
not affected by the very efficient cooling provided by the IR fine- 
structure lines of [O III] as in the high ionisation zones, where, in 
fact, the lowest temperatures are reached. The Z/Zq = 0.05 and 
0.2 models do not follow the same trend; here there is a monotonic 
shift of temperatures, whereby C models are hotter than the H and 
F models, regardless of ionic species zone. At intermediate metal- 
licities (Z/Zq = 0.4) we are in a transition case from the two sep- 
arate behaviours described above. At low metallicities, the cooling 
is dominated by collisional excitation of H Lya, which increases 
as one gets closer to the ionisation front, where the proportion of 
residual neutral hydrogen increases. This results in a outward de- 
creasing temperature, contrary to the high metallicity cases. 

Differences in the temperature structure of cases F and C can 
be understood as follows: at each point in the nebula the electron 
temperature is primarily determined by the average energy of the 
photons absorbed by H and He and by the cooling efficiency of the 
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Figure 2. H II region excitation sequence in terms of [N Il]/Ha and [S Il]/Ha vs.[0 IIl]/H/3 and R23. The grey points represent giant H II regions in galaxies 
taken from Gamett & Kennicutt (1994), Garnett et al. (1997), van Zee et al. (1998), Bresolin et al. (1999, 2004, 2005), and Kennicutt et al. (2003), spanning 
metallicities from 12 + log(0/H) ~ 8 to ~ 8.8., and metal-poor emission line galaxies from Izotov et al. (2006), which extend to 12 + log(0/H) ~ 7.2. The 
symbols and colors correspond to our models as defined in Table[2] Please refer to the on-line version of this paper for colour figures. 



ions. The latter is naturally less important at lower metallicities. 
However, for the intermediate to high metallicity cases, it is differ- 
ences in the distribution of the cooling ions rather than differences 
in the heating rates which explains the difference in temperature 
structure between case C and case H for higher and intermediate 
metallicities. 



3.2 Temperature Fluctuations 

We have calculated for the main ionic species in our models the 
formal values of the mean ionic temperatures, To, and of the tem- 
perature fluctuations, i 2 , according to the formalism of Peimbert 
(1967). Note that in the case of concentric sources, t 2 measures 
changes in radial gradients of the temperature rather than "tem- 
perature fluctuations" while in the cases of distributed sources it 
really represent a temperature fluctuation. For brevity we only list 
the values for H+, + and 2+ in TableO The errors quoted in 
the table are representative of the accuracy achieved by our mod- 
els. They contain contributions from the variance intrinsic to our 
Monte Carlo approach and the error introduced by using a finite 
grid to describe the ionised region. The errors were estimated by 
considering that for a fully spherically symmetric case, such as one 
with homogeneous gas density distribution and a central location 



for all ionising sources, the t 2 in an infinitely narrow spherical shell 
centred on the source of the ionising photons should be zero. We 
note that our errors increase with increasing metallicities, this is 
due to the larger temperature gradients occurring at higher metallic- 
ities causing the error contribution due to the finite grid description 
to increase. Clearly all errors could be reduced by increasing the 
number of energy packets and grid cells used in the simulations. 
However we note that our errors are always at least one order of 
magnitude smaller than the t 2 values and therefore do not affect 
our conclusions. Given the large number of 3D models run for this 
work (which vastly exceeds those finally presented here), we feel 
that a good balance between accuracy and computational expense 
was achieved. 

The value of t 2 (0 + ) is very hard to derive empirically. We 
know however of three H II regions where this measurement has 
been made: the Orion nebula (Esteban et al. 2004), M20 (Garcfa- 
Rojas et al., 2006) and M8 (Garcfa-Rojas et al. 2007). In those 
cases where t 2 (0 + ) cannot be derived empirically, it is generally 
assumed that t 2 (0 + ) = t 2 (0 2+ ); it is clear from the values listed 
in our table that this assumption is very often not verified and care 
should be taken to account for this in the error estimation from 
such studies. For our Z/Zq = 2 models, i 2 (0 2+ ) is always a fac- 
tor of 2 or more higher than t (0 + ), while for lower metallicities 
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Figure 3. H II region ionisation sequence parametrised in terms of oxygen and sulphur. The grey dots indicate the position of giant H II regions in spiral 
galaxies taken from Garnett & Kennicutt (1994), Garnett et al. (1997), van Zee et al. (1998), Bresolin et al. (1999, 2004, 2005), Kennicutt et al. (2003), 
spanning metallicities from 12 + log(0/H) ~ 8 to ~ 8.8. The top right panel also includes a few points from metal-poor emission line galaxies from Izotov et 
al. (2006), which extend to 12 + log(0/H) ~ 7.2; it was not possible to use this data for the other three panels, due to the lack of simultaneous detection of all 
the lines needed. The symbols and colors correspond to our models as defined in Table|2] Please refer to the on-line version of this paper for colour figures. 



t 2 (0 2+ ) becomes lower than t 2 (0 + ) sometimes by large factors 
(up to approximately 10). Finally, we note that the formal t 2 val- 
ues may diverge from the empirical ones (e.g. Kingdon & Ferland, 
1995; Zhang, Ercolano & Liu, 2007), that are generally based on 
the comparison of the electron temperature derived from the depth 
of the Balmer jump and the [O III] temperature, the values listed in 
Table |3TI are, nevertheless, sufficient to identify the effects of the 
stellar distribution on temperature fluctuations in the nebular gas. 

In Figure[5]we plot the t 2 values for H + and 2+ (dashed and 
solid lines, respectively) against metallicity for the spherical and 
shell models. Cases C, H and F are represented by black, red and 
green lines and symbols respectively. Unsurprisingly, the temper- 
ature fluctuations, which in this case are a direct consequence of 
large-scale temperature gradients, are larger for the high metallic- 
ity models. A full discussion of this effect and of the causes of the 
large scale temperature fluctuations in metal rich nebulae is given 
by Stasinska (1980) and Kingdon & Ferland (1995). What had not 
been noticed before, however, is that at very low metallicities, the 
t 2 values for H + and + rise again; this is due to the large tempera- 
ture gradient already shown in Figure|4] However, t 2 (0 2+ ) remains 
tiny, so that the overall effect on abundance determinations is ex- 
pected to be small. 

With regards to a comparison of the t 2 values obtained with 



the three different spatial distributions of sources (cases C, H and 
F), we first of all notice that for intermediate to high metallici- 
ties (Z/Zq > 0.4), case F models show a more isothermal gas 
(smaller t 2 ) than case H or C models. This is due to the fact that 
the contributions to the t 2 due to the true temperature fluctuations 
created by the two different stellar populations in case F (and in a 
smaller degree in case H) are completely washed out by the large- 
scale temperature gradients caused by the metals' cooling. This ef- 
fect vastly dominates at these metallicities, resulting in larger val- 
ues of t 2 being obtained by case C models, which, as discussed in 
the previous section (see Figure 4), have steeper temperature gradi- 
ents than their cases H and F counterparts. 



The above is further confirmed by the fact for very low metal- 
licities (Z/Zq < 0.2), case F models (green lines) show larger 
fluctuations than case C and H models. This is because at these 
low metallicities, case C, H and F all show similar large-scale tem- 
perature gradients (see Figure 4); the values of t 2 are thus larger 
for case F models where true temperature fluctuations are at play. 
However, once again we point out that no large effects are expected 
on abundance determinations due to the fact that t 2 (0 2+ ) is small 
in all cases. 
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Figure 4. Volume integrated electron temperatures weighted by ionic species. Top left: results for all models with Z/Zq = 2 (magenta); top right: results for 
all models at Z/Zq = 1 (black); middle left: results for all models at Z/Zq = 0.4 (red); middle right: results for all models at Z/Zq = 0.2 (green); bottom 
left: results for all models at Z/Zq = 0.05 (blue); bottom right: Results for all models at all metallicities. A legend for the symbols is given in Table[2] Please 
refer to the on-line version of this paper for colour figures. 



3.3 Ionic temperature relations 

We have shown that the temperature structure of models with cen- 
trally concentrated ionising sources, case C, may vary compared to 
those of similar models where the sources are distributed within the 
half volume (case H), which have partially overlapping Stromgren 
spheres, and to those with sources distributed within the full volume 
(case F), which have fully independent Stromgren spheres. Case 



C to case F variations in the temperature structure of the models 
may have implications for a number of ionic temperature relations. 
These scaling laws are often employed in abundance studies when 
observational data for a given ionic zone is missing. Some of the 
most popular relations we found in the literature include T e (S III) 
vs T e (Ar III), T e (0 III) vs T e (Ar III), T e (0 II) vs T e (N II). These 
ratios are expected to be around unity and we found them to be very 
little affected by the shifts in temperatures due to the spatial distri- 
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Figure 5. Mean ionic temperatures and temperature fluctuations as a function of gas metallicity for the homogeneous sphere and shell models. t 2 (0 2+ ) 
and t 2 ( H + ) are represented by solid and dashed lines, respectively. The lines are colour coded such that black, red and green are for cases C, H and F, 
respectively. Please refer to the on-line version of this paper for colour figures. 
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Figure 6. T c (0 II) versus T e (0 III) for Z/Z Q =1.0 and 2.0 left and Z/Z =0.05, 0.2 and 0.4 right models for cases C, H and F. The straight black solid line is 
the scaling law of Garnett (1992) and the blue dotted and red dashed curves are the relations given in Eqn 14 of Izotov et al. (2006) for 12 + log(0/H) = 8.2 
(red), 7.6 (blue dotted, thick) and 7.2 (blue dotted, thin). The symbols representing our models are described in Table|2] Please refer to the on-line version of 
this paper for colour figures. 
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Table 3. Mean temperatures and temperature fluctuations quantified in terms of T (H+), t 2 (H+), T (O+), t 2 (O+),T (O 2 +) and t 2 (0 2+ ) (Peimbert, 1967) 



model 


T (H+) 


t 2 {H+) 


T (O+) 


t 2 (o+) 


T (O+ 2 ) 


t 2 (0+ 2 ) 


no UTTft AC 

CSpnHO.05 


1 8240 


0.0178 


13350 


0.041 


18500 


0.007 


riipnrlU.Uj 


17690 


U.022o 


12830 


0.046 


18500 


0.007 


ropnrlU.UJ 


14850 


0.05 1 1 


1 1720 


0.035 


1 7 1 80 


A AAA 
U.UOV 






± 0.0006 




rb 0.003 




db 0.004 


L.iSpnrlU.2 


13140 


0.0034 


IxUoU 


0.0145 


13430 


0.0010 


ri.SpnriU.2 


13230 


0.0036 


11810 


0.0156 


13430 


0.0005 


roptlHU.2 


12100 


0.01 12 


1 1 160 


0.0141 


12860 


0.0010 






± U.UUU1 




i a f\r\f\A 
± U.UUU4 




1 n AAA1 

± U.UUU1 


L.iSpnrlU.4 


9850 


0.0063 


1 0620 


u.uuoy 


9790 


0.0053 


TIC UTTrt /I 

HSpnH0.4 


9860 


0.0038 


10500 


0.0060 


9790 


0.0029 


roptlHU.4 


9710 


0.0036 


9910 


0.0047 


9600 


A AAT A 
U.UU24 






I a aaac 

± 0.0005 




I A nnfn 
± 0.000/ 




I A AAAT 


Uopilrll.U 


6690 


0.024 


8130 


0.020 


6380 


0.012 


nopxiril.U 


6670 


0.021 


8110 


0.020 


o380 


0.009 


CCnliUl A 

ropnril.U 


7320 


0.018 


7990 


0.012 


o / 1U 


0.009 






± 0.006 




=b 0.002 




± 0.001 


Uopilrlz,U 


5000 


a i c./; 
O.Ijo 


61 10 


0.043 


2910 


0.130 


TIC Lin A 

HSpnHi.U 


4880 


0.152 


5980 


0.050 


2919 


0.096 


ropnH2.0 


5360 


0.075 


5890 


0.036 


3700 


0.084 






± 0.005 




rb 0.003 




1 0.003 


UonllrlU.Uj 


17510 


0.0129 


13770 


0.0437 


18110 


0.0041 


Ticuuun a c 

HinnHO.05 


17520 


0.0166 


13440 


0.0438 


18110 


0.0057 


ronnrlU.UD 


15200 


0.0369 


12410 


0.0348 


16840 


0.0U95 






I a nnnn 
± 0.0002 




I a Anni 
± 0.0002 




± 0.0001 


L.innrlU.2 


12940 


0.0033 


12390 


0.01 13 


13210 


0.0017 


riiSnnrlU.z 


13080 


0.0029 


12250 


0.0124 


13210 


0.0010 


ronilHU.2 


12280 


U.UUoZ 


1 1520 


0.0121 


12810 


0.0015 






± 0.0002 




± 0.0002 




± 0.0007 


CShhH0.4 


9780 


0.0091 


10570 


0.0089 


9710 


0.0078 


HShhH0.4 


9880 


0.0084 


10740 


0.0081 


9710 


0.0068 


FShhH0.4 


9810 


0.0060 


10130 


0.0064 


9620 


0.0046 






± 0.0003 




± 0.0003 




± 0.0001 


CShhHl.O 


6880 


0.030 


7980 


0.026 


6460 


0.017 


HShhHl.O 


6850 


0.029 


8010 


0.026 


6460 


0.016 


FShhHl.O 


7280 


0.021 


7940 


0.017 


6720 


0.012 






± 0.001 




±0.001 




± 0.001 


CShhH2.0 


5450 


0.086 


5980 


0.045 


3420 


0.09 


HShhH2.0 


5430 


0.090 


6000 


0.045 


3420 


0.09 


FShhH2.0 


5490 


0.070 


5950 


0.037 


3690 


0.08 






± 0.006 




± 0.006 




±0.01 



button of stars. This is not surprising, given that the ionic species 
involved are both 'high' (O III, S III, Ar III and N III) or Tow' (O II 
and N II), and therefore while the absolute temperature values in 
each case may shift to higher or lower values, the resulting ratios 
remain virtually unaffected. 

Scaling laws have also been derived for T e (0 III) vs T e (S III) 
and T e (0 III) vs T e (0 II), by Garnett (1992), based on his own 
grids of photoionisation models and those by Stasinska (1982). 
More recently, Izotov et al. (2006) produced somewhat different 
relations between T e (0 II) and T e (0 III), and between T e (S III) 
and T e (0 II) (their Eqs. 14 and 15, valid only for metal poor cases 
to 12 + log(0/H) ~ 8.2), based on a set of up-to-date (but still 
spherically symmetric) photoionisation models that reproduce the 
observed trends of metal-poor galaxies. They found that the ob- 
served values of T e (0 II), T e (0 III) and T e (S III) in a large sample 
of H II galaxies reproduce the theoretical relations, but with a large 
scatter (not only attributable to observational errors in the case of 
T e (0 II), T e (0 III), see their Figures 4a and b). 

Figure [6] shows the behaviour of our Z/Z© =1.0 and 2.0 (left 



panel, black and magenta symbols respectively) and Z/Z© =0.05, 
0.2 and 0.4 (right panel, blue, green and red symbols, respectively) 
metallicity models for cases C, H and F with regards to the T e (0 II) 
versus T e (0 III) scaling laws of Garnett (1992, black solid line in 
both panels) and Izotov et al. (2006, coloured curves in the right 
panels). The symbols are as described in Table[2] Our lowest metal- 
licity bin Z/Z = 0.05 (corresponding to 12 + log(0/H)= 7.39) falls 
in between Izotov's 12 + log(0/H)= 7.2 and 7.6 metallicity bins. 
Their T(0 III) versus T(0 II) relations given in their Eqns 14 are 
therefore plotted in blue in the right panel of our Figure [6] with 
the thinner and thicker curves indicating the 7.2 and 7.6 metal- 
licity bin, respectively. Our Z/Z© = 0.4 metallicity bin (red sym- 
bols) corresponds to 12 + log(0/H)= 8.29, slightly above Izotov's 
12 + log(0/H) = 8.2 bin which, nevertheless is represented by the 
red line in Figure[6] 

A 'mixed' relation such as T e (0 III) vs T e (0 II) is much more 
affected by the fact that low and high ionic temperatures shift in 
opposite directions for intermediate (red points; Z/Z© = 0.4) and 
high metallicities (magenta and black points; Z/Z© = 2 and 1). This 
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Table 4. Oxygen abundances derived applying various metallicity indicators to the integrated emission line spectra calculated for our models (See text for 
calibrations used). All abundances are given in logarithmic scale, 12 + log(0/H). The input values of 12 + log(0/H) (i.e. 'the right answer') are 8.99 (Z/Zq = 2), 
8.69 (Z/Zq = 1), 8.29 (Z/Z = 0.4), 7.99 (Z/Z = 0.2) and 7.39 (Z/Z© = 0.05). A |, the largest difference between the abundance derived from C, H and F 
cases are given for each model trio. The value averaged over all models, |< A >|, is given in the last row of the table. 



model 


023 


O3N2 


N 2 


S23 


S3O3 


Ar 3 3 


model 


R23 


O3N2 


N 2 


S23 


S3O3 


Ar 3 3 


CSp2.0 


8.67 


8.78 


8.44 


7.92 


8.77 


8.83 


CSh2.0 


8.62 


8.77 


8.50 


8.09 


8.75 


8.82 


HSp2.0 


8.69 


8.78 


8.43 


7.87 


8.76 


8.81 


HSh2.() 


8.62 


8.77 


8.50 


8.08 


8.75 


8.82 


FSp2.0 


8.67 


8.73 


8.52 


8.17 


8.67 


8.78 


FSh2.0 


8.63 


8.76 


8.52 


8.16 


8.72 


8.81 


A 


0.02 


0.05 


0.09 


0.30 


0.10 


0.05 


1 A | 


0.01 


0.01 


0.02 


0.08 


0.03 


0.01 


CSpl.O 


8.74 


8.27 


8.17 


7.76 


8.04 


8.48 


CShl.O 


O.DO 


8.30 


8.24 


"7 CiA 


5. ID 


6.55 


HSpl.O 


8.75 


8.27 


8.16 


7.73 


8.02 


8.46 


HShl.O 


8.67 


8.30 


8.23 


7.92 


8.14 


8.54 


FSpl.O 


8.49 


8.40 


8.43 


8.34 


8.24 


8.62 


FShl.O 


8.51 


8.38 


8.39 


8.23 


8.24 


8.61 


A 


0.26 


0.13 


0.27 


0.61 


0.22 


0.16 


1 A | 


0.16 


0.08 


0.16 


0.31 


0.10 


0.07 


CSp0.4 


8.52 


8.07 


7.87 


7.93 


8.28 


7.96 


CSh0.4 


8.52 


8.09 


7.97 


8.49 


8.10 


7.77 


HSp0.4 


8.52 


8.06 


7.85 


7.92 


8.22 


7.91 


HSh0.4 


8.50 


8.10 


7.98 


8.45 


8.08 


7.73 


FSp0.4 


8.43 


8.24 


8.15 


8.22 


9.02 


8.25 


FSh0.4 


8.41 


8.23 


8.19 


8.92 


8.27 


8.12 


A 


0.09 


0.18 


0.30 


0.30 


0.80 


0.34 


|A| 


0.11 


0.14 


0.22 


0.47 


0.19 


0.39 


CSp0.2 


7.86 


7.97 


7.77 


7.99 


7.67 


7.67 


CSh0.2 


7.84 


7.99 


7.79 


8.18 


7.87 


7.92 


HSp0.2 


7.87 


7.97 


7.77 


8.01 


7.66 


7.66 


HSh0.2 


7.86 


7.97 


7.77 


8.10 


7.79 


7.82 


FSp0.2 


7.68 


8.15 


8.05 


8.68 


8.10 


8.26 


FSh0.2 


7.72 


8.13 


8.02 


8.61 


8.11 


8.26 


A 


0.19 


0.18 


0.28 


0.69 


0.44 


0.60 


1 A| 


0.14 


0.16 


0.25 


0.51 


0.32 


0.44 


CSp0.05 


7.48 


7.92 


7.51 


7.29 


7.39 


6.30 


CSh0.05 


7.44 


7.93 


7.51 


7.38 


7.58 


6.67 


HSp0.05 


7.44 


7.95 


7.56 


7.37 


7.49 


6.51 


HSh0.05 


7.44 


7.94 


7.53 


7.38 


7.57 


6.63 


FSp0.05 


7.15 


8.12 


7.77 


7.73 


7.91 


7.48 


FSh0.05 


7.21 


8.09 


7.74 


7.69 


7.93 


7.49 


A 


0.33 


0.20 


0.26 


0.44 


0.52 


1.18 


|A| 


0.23 


0.16 


0.23 


0.31 


0.36 


0.86 


<| A[> 


0.18 


0.15 


0.24 


0.47 


0.41 


0.47 


<l A |> 


0.12 


0.11 


0.18 


0.32 


0.23 


0.45 



will contribute to the scatter noticed observationally for this rela- 
tion. For example, Kennicutt, Bresolin & Garnett (2003) presented 
a study of those H II regions in M 101 for which direct measure- 
ments of the nebular auroral lines could be obtained. They found 
the T e (0 III) vs T e (S III) relation to be matched closely by their ob- 
servations, whereas the T e (0 III) and T^O II) temperatures turned 
out to be rather uncorrelated. It is also worth noting at this point tht 
another factor that could contribute to the scatter in the T e (0 III) 
vs T e (0 II) relation is the dependency of T e (0 II) on the electron 
density (Perez-Montero & Diaz 2003). The scatter found by our 
models is for a given electron density, and therefore it may be even 
larger for a different sample covering a wider range of densities. 
The ionising flux distribution is certainly a factor that contributes 
to weakening the correlation, however, as pointed out by Stasinska 
(2005), at high metallicities, the temperature derived from the [O II] 
line ratio may by strongly in error due to the contribution of recom- 
bination from ++ . 



3.4 Abundance diagnostics 

We have assessed that the geometrical distribution of stars within 
an H II region plays a role in the temperature structure of the gas, 
the magnitude of these effects will clearly be dependent on the to- 
tal mass of the stars and the concentration level. It is important 
now to verify the robustness of commonly used metallicity indica- 
tor against these temperature shifts. In particular we are interested 
in identifying systematic trends with stellar distribution rather than 
absolute errors. Table [4] lists the empirical oxygen abundances ob- 
tained using 6 different indicators. We used the following calibra- 
tions for the 6 indicators: (1) O23 - Pilyugin (2000, 2001b), (2) 
3 N 2 - Stasinska (2006), (3) N 2 - Pettini & Pagel (2004), (4) S23 



- Perez-Montero & Diaz (2005), (5) S3O3 - Stasinska (2006), (6) 
Ar3C>3 - Stasinska (2006). It is worth reiterating at this point that 
metallicity indicators should, and generally are, only used for sta- 
tistical studies as the error on a single region may be very large; the 
models presented in this work do not attempt to cover the full pa- 
rameter space occupied by H II regions and therefore do not aim at 
identifying the most accurate indicator in absolute terms. Here we 
are mainly interested in studying a possible systematic error in the 
derived abundances introduced by the 3D stellar distribution. For 
this reason for each trio of models (C, H and F for a given density 
distribution and metallicity) we compute A, the largest difference 
between the oxygen abundance in log units derived from cases C, 
H and F. 

The results are summarised in Table [4] where the mean val- 
ues < A > are also listed. For all metallicity indicators apart from 
O23, there is a clear trend for higher metallicity being derived from 
models ionised by fully distributed sources, case F, than from case 
C and H models. The reverse is true for O23, where metallicities de- 
rived from case F models are smaller than those derived from case 
C and H models. The A values at a given metallicity vary from 
one indicator to the next, however they are rarely below ~0. 1 dex 
(only for the highest metallicity case), with more representative val- 
ues around 0.3 dex for the spherical case, but often larger than 0.4 
(and larger than 1.1 dex in one case). These deviations are slightly 
smaller for the shell density distributions; this is obvious as even 
for F cases most of the ionising radiation in these models will be 
emitted from the central cavity, reducing the differences between 
C and F cases. We note that in some cases the values given by the 
case H models are lower than those given by the respective case C 
models (rather than being equal or in the middle between C and F). 
We have analysed these deviations statistically and can confirm that 



Temperature in H II regions 1 1 



the small differences are simply due to the variance of our Monte 
Carlo models and do not bear any physical significance. 

The reason for the systematic effect we see in the empirical 
abundance determinations can be understood by re-examining Fig- 
ure [3] where the ionisation sequence of our model is parametrised 
in terms of oxygen and sulphur emission lines. The C (small sym- 
bols), H (medium symbols) and F (large symbols) form a sequence, 
with the case F models generally showing a lower ionisation param- 
eter than case C and H models. The differences in the abundances 
derived from case C to F with the various indicators clearly reflect 
the different dependence of each indicator on the ionisation param- 
eter. This can be simply shown for an idealised system. The ioni- 
sation parameter of a pure-hydrogen spherical volume of gas with 
number density Nh ionised by a source emitting Qh° hydrogen- 
ionising photons per second is defined as U = Q h o/(4tyNhR 2 ), 
where R s is the Stromgren radius. In such a system, from the ion- 
isation balance equation (e.g. Osterbrock, 1989, eqn 2.19) R s is 
directly proportional to Q^o an d to N^ 2 ^ , neglecting the tem- 
perature dependence of the hydrogen recombination coefficient. 
For simplicity we compare the ionisation parameter of the system 
above, U c , ionised by only one source, to that of a similar sys- 
tem (F) ionised by N s identical sources each emitting Q H olN s 
hydrogen-ionising photons per second. With these assumptions it 
is easy to show that the ionisation parameter, U F ', for system F, 
measured at the Stromgren radius of each individual source is sim- 
ply related to U c by 

U F = U c -N^ (2) 

and therefore always smaller than U c for N s > 1 

The magnitude of the errors in the metallicity derived by the 
strong line methods are significant. It is, however, true that the ef- 
fects reported here represent a worst case scenario, and our extreme 
assumptions on stellar populations create a large dispersion in the 
resulting ionisation parameters, which in some cases exceeds the 
observed range, as seen in Figures [2] and [5] Aside from the mag- 
nitude of the A values, however, a worrying aspect is the fact that 
the discrepancies between cases C, H and F are systematic. This 
can have an impact on galactic metallicity gradients determined 
via strong line methods, if these are calibrated via photoionisation 
models. In fact, if compact clusters (close to case C) and loose as- 
sociations (close to case F) are randomly distributed throughout a 
given galaxy, then the systematic errors would only cause a larger 
scatter in the observed metallicities, but, given sufficient number 
statistics, they would not affect the measured metallicity gradient. 
However if the ratios of compact clusters to loose associations is 
somewhat dependent on the galactocentric distance, then the sys- 
tematic errors due to stellar geometrical distributions may indeed 
introduce a bias on the measured galactic metallicity gradient, if 
the abundances are obtained from strong line methods calibrated 
on ab-initio models which do not reproduce the observed exci- 
tation of H II regions. For example, recent work by Rosolowsky 
et al. (2007), presenting high resolution molecular gas maps of 
M33, showed a truncation in the mass distribution of giant molec- 
ular clouds (GMCs) at a galactocentric distance of 4 kpc. A recent 
study on the demographics of young star-forming clusters in M33 
by Bastian et al. (2007) also shows the same cut-off at 4 kpc for the 
clusters detected. We could interpret this as tentative evidence of 
different star formation environments from the centre to the edge 
of M33, however, we prefer to postpone this discussion until more 
compelling observational evidence becomes available. 



4 CONCLUSIONS 

Following our theoretical investigation on the effects of the spatial 
configuration of ionisation sources on the temperature structure of 
H II regions, we summarise our conclusions as follows: 

(i) For intermediate to high metallicities (0.4 < Z/Zq < 2), 
for a given gas density distribution, abundance and ionising spectral 
shape and intensity, a model with a central concentration of stars 
(case C) will result in higher ionic temperatures for high ionisation 
species (0 2+ , S 2+ etc.), compared to the same model with stars 
fully distributed within the volume (case F). The opposite is true for 
Tow' ionisation species (e.g. + , N + ). This results in a shallower 
gradients in the electron temperature distribution across the ionic 
species zones. 

(ii) Low metallicity models (Z/Zq < 0.2) do not show the tem- 
perature inversion from low to high ionic species zones, rather a 
shift in the temperature is experienced by all ionic species zones, 
resulting in case F models being cooler than case C and H models. 

(iii) At intermediate to high metallicity (Z/Zq > 0.4), mod- 
els with stars distributed within the full volume are more isothermal 
(show lower t 2 values) than the same models with a central concen- 
tration of stars. These temperature "fluctuations" obtained for case 
C models are a simple consequence of a large temperature gradi- 
ent. Multiple ionising sources of different temperatures at central or 
non-central locations do not produce significant temperature fluc- 
tuations in the ionised gas of models with Z/Zq > 0.4. 

(iv) At low metallicities (Z/Zq < 0.2), models with stars 
distributed within the full volume (case F) show larger t 2 values 
than the same models with a central concentration of stars. Here 
we are seeing the effects of true temperature "fluctuations" for case 
F models. The magnitude of t 2 (0 2+ ) remains however too small 
to have any significant effect on derived abundances. 

(v) Multiple ionising sources of different temperatures at cen- 
tral or non-central locations are not the cause of significant 
temperature fluctuations in the ionised gas of our models with 
0.05 < Z/Zq < 2. 

(vi) The relation t 2 (0 + ) = t 2 (0 2+ ), often used in empirical 
studies, is NOT verified by our models. Extreme care should be 
taken to account for the uncertainties introduced by the use of 
this relation in studies seeking to apply corrections to CEL-derived 
abundances making use of an empirical estimation of temperature 
fluctuations. For our Z/Zq = 2 models, t 2 (0 2+ ) is always a fac- 
tor of 2 or more higher than t 2 (0 + ), while for lower metallicities 
t 2 (0 2+ ) becomes lower than t 2 (0 + ) sometimes by large factors 
(up to approximately 10). 

(vii) For intermediate to high metallicity models, electron tem- 
peratures in the 2+ and + ionisation zones are shifted in oppo- 
site directions, contributing to the scatter observed in the T e (0 ) 
versus T e (0 + ) relation. We confirm that H II region abundances 
derived on the basis of the T e (0 II) alone should, therefore, be con- 
sidered highly uncertain. 

(viii) Metallicity indicators calibrated by grids of spherically 
symmetric photoionisation models may suffer a systematic bias, 
due to their dependence on the ionisation parameter of the system. 
For the same input parameters case F models will always result in 
smaller ionisation parameters than case C and H models. The er- 
rors estimated in this work (typically 0.3 dex, but larger in some 
cases) are likely to represent the worst-case scenario, but never- 
theless their magnitude and their systematic nature does not allow 
them to be ignored. 
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